En esta unidad se abordan los principios de los modelos de regresión aplicados a las series temporales. Se introducen los conceptos de estacionariedad, autocorrelación y autocovarianza, y cómo estos son fundamentales para trabajar con modelos de series temporales.
Regresión simple y múltiple en el contexto de las series temporales
Estacionariedad en series temporales
Funciones de autocorrelación y autocovarianza
Análisis de la autocorrelación y parcial autocorrelación en series temporales
Modelos autoregresivos (AR) y medias móviles (MA)
Regresión lineal en series de tiempo
Concepto clave
Podemos aplicar el concepto de regresión para pronosticar una variable \(y_t\) asumiendo que tiene una relación esperada con otra en el tiempo \(x_t\).
import numpy as npimport matplotlib.pyplot as pltimport statsmodels.api as smx = np.random.normal(50, 10, 100)y = np.random.normal(2, 10, 100) +2* x + np.random.normal(0, 1, 100)X = sm.add_constant(x)lm = sm.OLS(y, X).fit()y_hat = lm.predict()plt.style.use('ggplot')fig, ax = plt.subplots()ax.scatter(y, x, color ="gray")ax.plot(y_hat, x, color ="red")plt.title('Regresión lineal')plt.xlabel('X')plt.ylabel('Y')fig.show()
Estimador Minimos Cuadrados Ordinarios (MCO)
Definición
El estimador MCO permite buscar los coeficientes \(\beta_k\) de una regresión tal que minimiza la suma de los errores (el valor más pequeño) de la regresión.
Donde las dimensiones del sistema \(n\) filas, \(k\) columnas : \(\mathbf{y} = n \times 1\), \(\mathbf{X} = n \times k\), \(\mathbf{\beta} = k \times 1\) y \(\mathbf{\varepsilon} = n \times 1\)
Teniendo en cuenta que MCO minimiza el error cuadrático en términos matriciales \(\sum_{t=1}^T \varepsilon_t^2 = \mathbf{\varepsilon}'\mathbf{\varepsilon}\)
Finalmente, tenemos que el estimador MCO es :\[
\mathbf{\hat{\beta}} = \mathbf{X}'\mathbf{X} ^ {-1} \mathbf{y}'\mathbf{X}
\]
Nota
Si tenemos un sistema lineal de la forma \(\mathbf{A}\mathbf{x} = \mathbf{b}\). La solución del sistema no homogéneo es \(\mathbf{x} = \mathbf{A}^{-1} \mathbf{b}\) sí y solo sí \(\mathbf{A}\) es una matriz invertible no singular.
Nota
Propiedad la transposición de matrices:
\((\mathbf{ABC})' = \mathbf{C'B'A'}\).
Si \(\mathbf{A}\) es simétrica entonces \(\mathbf{A}'\) = \(\mathbf{A}\)
Si \(\mathbf{A}\) es inversible sí y solo sí \(\mathbf{A}'\mathbf{A} = \mathbf{I}\)
Derivada de una matriz \(a'x\), \(\frac{\partial (a'x)}{\partial x} = a\)
Derivada de una matriz cuadrática, \(x'Ax\)\(\frac{\partial (x'Ax)}{\partial x} = 2Ax\)
La interpretación de los \(\beta\)
Si tanto \(Y\) como \(X\) son continuos \[\beta_k = \frac{\Delta Y}{\Delta X_k}\]\(\beta_k\) = El cambio promedio esperado en \(Y\) ante un cambio en \(X\)
Si \(Y\) es númerica y \(X\) es dicotómica \[\beta_k = E(Y|X = 1) - E(Y|X = 0)\]\(\beta_K\) = Es la diferencia promedio esperada en \(Y\) cuando \(X = 1\) respecto a \(X = 0\)
Supuestos del modelo lineal
Las observaciones provienen de una muestra aleatoria de una población
Independecia : \(Y_i, .. Y_n\) son independientes
Linealidad : El modelo de regresión es lineal en los parámetros, aunque puede o no ser lineal en las variables.
Lineal \[Y = \beta_0 + \beta_1 \log X\]
No lineal \[Y = \beta_0 + 1/\beta_1 X\]
Homocedasticidad : El valor medio del error es igual a cero y La varianza del término de error, o de perturbación, es la misma sin importar el valor de X. \[
\varepsilon_i \sim N(0, \sigma^2_{\varepsilon})
\]
No hay autocorrelación entre las perturbaciones \[
Cov(\varepsilon_i, \varepsilon_j) = 0
\]
Exogeneidad : Valores de X independientes del término de error \(Cov(\varepsilon_i, \mathbf{X}) = 0\)
Colinealiad perfecta : no existe correlación perfecta entre las regresoras incluidas en el modelo de regresión \(X_1 = 2X_2\)
Rango completo : El número de observaciones n debe ser mayor que el número de parámetros por estimar \(N > X_k\)
Propiedades del estimador MCO
El supuesto de \(\varepsilon \sim N(0, \sigma^2_{\varepsilon})\) tiene las siguientes implicaciones para los parámetros estimados por MCO
Son insesgados\[
E(\hat{\beta_k}) = \beta_k
\]
Tienen varianza mínima. En combinación con 1, esto significa que son estimadores insesgados con varianza mínima, o eficientes. \[
Var(\hat{\beta_0}) = \frac{\sum X^2_i}{n\sum(X_i - \overline{X})^2}, Var(\hat{\beta_1}) = \frac{\sigma^2}{\sum(X_i - \overline{X})^2}
\]
Presentan consistencia; es decir, a medida que el tamaño de la muestra aumenta indefinidamente, los estimadores convergen hacia sus verdaderos valores.
Ejercicio 1. Gráficar la convergencia de los parámetros estimados por MCO a sus verdaderos valores.
Supongamos que los valores de \(y\) siguen esta función lineal simple \[
y = 12 + 20 x + \varepsilon
\]
Donde \(\beta_0 = 12\), \(\beta_1 = 20\), \(X \sim N(12, 24)\) y \(\varepsilon \sim N(0, 1)\)
Simular la convergencia de los parámetros estimados por MCO según incrementa la muestra.
Code
import numpy as npimport matplotlib.pyplot as pltimport statsmodels.api as smn_sim = np.linspace(10, 10000, 100).astype(int)true_beta = {"b0" : 12, "b1" : 20}x_stats = {"mu" : 12, "sigma" : 24}beta_hat = {"b0" : [], "b1" : []}# Generar simulacionfor n in n_sim: np.random.seed(42+ n) x = np.random.normal(x_stats['mu'], x_stats['sigma'], n) y = true_beta['b0'] + true_beta['b1'] * x + np.random.normal(0, 1, n) X = sm.add_constant(x) lm = sm.OLS(y, X).fit() beta_hat['b0'].append(lm.params[0]) beta_hat['b1'].append(lm.params[1])plt.style.use('ggplot')fig, axes = plt.subplots(1, 2, figsize=(12, 5))colors = ['#d95f02', '#1b9e77']# Iteramos sobre los resultados para graficarfor i, (key, value) inenumerate(beta_hat.items()): axes[i].plot(n_sim, value, label=f'Estimado {key}', color=colors[i], marker='o') axes[i].axhline(true_beta[key], color='black', linestyle='--', alpha=0.6, label=f'Real {key}') axes[i].set_title(f'Convergencia de {key}') axes[i].set_xlabel('Tamaño de muestra (n)') axes[i].set_ylabel('Valor del parámetro') axes[i].legend()plt.tight_layout()plt.show()
Medición del ajuste global del modelo lineal
La bondad del ajuste mide cuán “bien” se ajusta la línea de regresión a los datos. Se calcula en el coeficiente de determinación \(R^2\).
Esta mide la proporción o el porcentaje de la variación total en \(Y\) explicada por el modelo de regresión. \[
R^2 = \frac{\sum(\hat{y}_{t} - \bar{y})^2}{\sum(y_{t}-\bar{y})^2}
\]
Sus valores \(R^2 > 0\) en el rango \(0 \leq R \leq 1\)
Note
Cuando la regresión es simple, es decir con una sola variable regresora \(x_t\) el coeficiente de determinación \(R^2\) es el mismo que la correlación lineal. Recordemos que, la correlación lineal mide la fuerza y el sentido de una relación entre dos variables numéricas. \[
r = \frac{\sum (x_{t} - \bar{x})(y_{t}-\bar{y})}{\sqrt{\sum(x_{t}-\bar{x})^2}\sqrt{\sum(y_{t}-\bar{y})^2}}
\]
Significancia de los parámetros
Los parámetros estimados por MCO se le verifica su significancia a través de una prueba de hipótesis
\[
H_0: \beta_j = 0
\]\[
H_a: \beta_j \neq 0
\]
\[
t_{1 - \alpha/2; (n - k)} = \frac{\beta_j - \hat{\beta_j}}{\sqrt{Var(\hat{\beta_j})}}
\] - Regla de descisión : Si \(P(t_{1 - \alpha/2; (n - k)}) < 0.05\) consideramos que hay evidencia suficiente contra \(H_0\), por lo que se acepta \(H_a\). \[
\text{IC95 \%} = \hat{\beta_j} \pm t_{1 - \alpha/2; (n - k)} \sqrt{Var(\hat{\beta_j})}
\]
Ejemplo 2 : El modelo macroeconómico del crecimiento del consumo en Estados Unidos
Queremos modelar cómo cambia el consumo real (compra de bienes y servicios) de una economía según variables macroeconómicas importantes. El consumo es uno de los componentes más importantes de una economía.
# | eval: false# Cargar librerías necesariaslibrary(tidyverse)library(fable)library(tsibble) # 'tsibble' para series de tiempo# 1. Cargar y preparar datosus_change <- read_csv("datos/US_change.csv") %>% mutate(ds = yearquarter(ds)) %>% as_tsibble(index = ds)# 2. Visualizaciónggplot(us_change, aes(x = ds)) + geom_line(aes(y = Income, color = "Income")) + geom_line(aes(y = y, color = "Consumption")) + labs(x = "Quarter [1Q]", y = "% change", color = "Variable") + theme_minimal()# 3. Modelo de Regresión Lineal (MCO)# Usamos TSLM (Time Series Linear Model) de la librería 'fable'fit <- us_change %>% model( mco = TSLM(y ~ Income + Production + Savings + Unemployment) )# 4. Resumen del modeloreport(fit)
Evaluando el modelo MCO
Tras seleccionar las variables de regresión y ajustar el modelo, es necesario graficar los residuos para verificar que se han cumplido los supuestos del modelo.
Recordemos que los errores de estimación del modelo : \[
\varepsilon_t = y_t - \hat{y}_t \\
\varepsilon_t = y_t - \hat\beta_{0} - \hat\beta_{1} x_{1,t} - \hat\beta_{2} x_{2,t} - \cdots - \hat\beta_{k} x_{k,t} \\
\sum_{t=1}^{T}{e_t}=0 \quad\text{and}\quad \sum_{t=1}^{T}{x_{k,t}e_t}=0\qquad\text{for all $k$}. \\
\text{Para todo} t=1,\dots,T
\]
Linealidad de los paramentros
La relación de las variables regresoras deben ser lineales con la variables explicada. Gráficamente, no debemos ver un patron funcional entre los errores y los valores predichos del modelo.
Distribución de los errores son normales
Al hacer un histograma de los errores del modelo buscamos observar distribuciones sesgadas o no propia de la distibución normal.
Para comprobar formalmente podemos ayudarnos de un test de normalidad evalúa si un conjunto de datos puede considerarse proveniente de una distribución normal. Usualmente se utiliza Shapiro–Wilk.
Residuos del modelo con varianza constante (homocedasticos)
Los errores del modelo deben distibuirce de forma constante. Gráficamente, no debemos ver un patron entre los errores y los valores predichos del modelo.
Residuos no está correlacionados
Es usual que los datos temporales tenga correlación con datos anteriores debido a que sus valores se repiten en el tiempo, como el componente estacional. En este caso, el supuesto de no autocorrelación de los errores del modelo MCO se viola y los estimadores son ineficientes pero aún son insesgados.
Es un gráfico que mide la correlación que existe entre el valor de una serie de tiempo y su rezago o “lag” \(K\)\[
r_{k} = \frac{\sum\limits_{t=k+1}^T (y_{t}-\bar{y})(y_{t-k}-\bar{y})}
{\sum\limits_{t=1}^T (y_{t}-\bar{y})^2},
\]
Cuando los datos presentan una tendencia, las autocorrelaciones para los rezagos (lags) pequeños tienden a ser grandes y positivas, debido a que las observaciones actuales y anteriores son cercanos entre sí. Por lo tanto, la función de autocorrelación (ACF) de una serie con tendencia tiende a mostrar valores positivos que decrecen lentamente a medida que los desfases aumentan.
En cambio, cuando los datos son estacionales, las autocorrelaciones serán mayores para los rezagos estacionale (semanales, trimestrales, anuales, etc) que para otros rezagos.
¿Puedes análizar estos ACF? ¿Tienen tendencia y estacionalidad?
Test de Durbin-Watson para la autocorrelación serial
Es una prueba estadística de prueba que se utiliza para detectar la presencia de autocorrelación en los residuos de un análisis de la regresión.
\[
d = \frac{\sum_{t=2}^{T} (e_t - e_{t-1})^2}{\sum_{t=1}^{T} e_t^2} \\
0 \leq d \leq 4
\]
Regla de desición
d = 2 indica que no hay autocorrelación
d < 2 indica correlación serial positiva
d > 2 indica correlación serial negativa
Test de Ljung-Box para la autocorrelación
Además de observar el gráfico de la ACF, también podemos realizar una prueba más formal de autocorrelación considerando un conjunto completo de rezagados \(r_k\) como un grupo, en lugar de tratar a cada uno por separado.
Important
Ruido Blanco (White Noise): Una serie donde los datos son totalmente independientes e idénticamente distribuidos. En la ACF, esto se ve como picos que no superan los límites de confianza.
H0 : Los residuos no están correlacionados hasta el \(\ell\) resagos (no son distintos de un proceso ruido blanco)
Ha : Al menos uno de los residuos está correlacionado hasta el \(\ell\) resago
Regla de desición : Si \(P(Q^*) > 0.05\) la serie sigue un proceso ruido blanco.
Soluciones al problema de autocorrelación
Para una correlación serial positiva, considere agregar rezagos de la variable dependiente y / o independiente al modelo. \[
y_t = \beta_0 + \beta_1 x_t + \beta_2 y_{t-1} + \varepsilon
\]
Para una correlación serial negativa, verifique que ninguna de sus variables esté sobrediferenciada. \[
\Delta \Delta x_t
\]
Para la correlación estacional, considere agregar variables ficticias estacionales al modelo. \[
y_{t} = \beta_{0} + \beta_{1} x_t + \beta_{2}d_{2,t} + \beta_3 d_{3,t} + \beta_4 d_{4,t} + \varepsilon_{t}
\]
Ejercicio 3 :
Hacer el diagnostico del modelo macroeconómico del crecimiento del consumo en Estados Unidos. ¿Se viola algún supuesto del MCO?
Solución
Code
from statsmodels.graphics.tsaplots import plot_acffrom statsmodels.stats.diagnostic import acorr_ljungboxfrom statsmodels.stats.stattools import durbin_watsondef plot_diagnostics(modelo, n_lags=None, target_col="y", time_axis=None):ifnot"statsmodels"instr(type(modelo)):raiseTypeError("El objeto debe ser un modelo ajustado de statsmodels.")# Extraemos residuos directamente del objeto de resultados resid = modelo.resid# Si no se provee un eje de tiempo, usamos el índice de los residuosif time_axis isNone: time_axis = resid.index fig, axes = plt.subplot_mosaic( [["resid", "resid"], ["acf", "hist"]], figsize=(12, 8), constrained_layout=True )# Gráfico de Residuos ax = axes["resid"] ax.plot(time_axis, resid, marker='.', linestyle='-', alpha=0.7) ax.axhline(0, color='black', linestyle='--', linewidth=1) ax.set(title=f"Innovation Residuals ({target_col})", ylabel="Residual")# Gráfico ACF ax = axes["acf"] plot_acf(resid, ax=ax, zero=False, lags=n_lags, auto_ylims=True) ax.set(title="ACF Plot", xlabel="Lag", ylabel='ACF')# Histograma ax = axes["hist"] ax.hist(resid, bins=20, edgecolor='white', color='steelblue') ax.set(title="Histogram", xlabel="Residual value", ylabel="Count") plt.show()# Ljung-Box (Prueba global) lb_test = acorr_ljungbox(resid, lags=n_lags, return_df=True) min_p_value = lb_test['lb_pvalue'].min() is_white_noise = min_p_value >0.05# Durbin-Watson (Autocorrelación de primer orden) dw_stat = durbin_watson(resid) result_test = (f"PRUEBA DE DIAGNÓSTICO DE RESIDUOS\n"f"{'-'*30}\n"f"Durbin-Watson: {dw_stat:.2f} (Objetivo: ≈ 2.0)\n"f"Ljung-Box Test (up to lag {n_lags}):\n"f" - Min p-value: {min_p_value:.4f}\n"f" - Conclusión: {'Ruido Blanco'if is_white_noise else'Autocorrelación'}" ) print(result_test)# Ejecución de la funciónplot_diagnostics(ml, n_lags=20, target_col="Consumption", time_axis=us_change.index)
PRUEBA DE DIAGNÓSTICO DE RESIDUOS
------------------------------
Durbin-Watson: 2.22 (Objetivo: ≈ 2.0)
Ljung-Box Test (up to lag 20):
- Min p-value: 0.0211
- Conclusión: Autocorrelación
Predictores temporales
Es usual que usemos caracteristicas de la serie temporal en la regresión tales como :
Tendencia
La tendencia \(T\) se puede modelar incluyendola como variable explicativa en el modelo \[
y_{t}= \beta_0+\beta_1t+\varepsilon_t
\]
Variables dummy
Las variables dummy que son variables que solo toman valores \({0, 1}\) se utilizan para :
Crear variables estacionales en los datos, tal que capturen patrones especificos de esos momentos. Puedes capturar picos específicos (como Navidad o vacaciones de verano) que una tendencia lineal simple ignoraría. \[
y_{t} = \beta_{0} + \beta_{1} t + \beta_{2}d_{2,t} + \beta_3 d_{3,t} + \beta_4 d_{4,t} + \varepsilon_{t}
\]
Caution
¡Importante! : Debes dejar siempre una categoría de comparación por lo que la cantidad de dummys a incluir será \(d-1\) variables. En pandas usamos get_dummies(drop_first = True). Si omites este paso, existe colinealidad perfecta y el modelo no se puede estimar.
Variables dummy estacionales
Son variables dummy pero que marcan eventos estacionales en la serie, tales como días de la semana, meses de un año, trimestres, etc. Es importante, dejar una variable por fuera del modelo por el riesgo de colinealidad perfecta.
Un problema que podemos enfrentar es cuando la serie tiene una estacional alta como las series semanales, ya que por año son 52 semanas, lo cual termina en 51 variables dummy estacionales. Recordemos que para cada variable el estimador MCO requiere estimar un \(\hat{\beta}\) lo cual es más dificil en la medida que incluimos más variables con los mismos datos.
Series de Fourier
Es una alternativa a las dummy estacionales, en especial en series con estacionalidad alta. Jean-Baptiste Fourier un matemático Fránces describió que las funciones seno y coseno pueden usarse para aproximarse a cualquier serie periodica. Aquí vemos las primeras pares de la serie de seno y cosenos para una estacional de una semana :
El máximo de pares de seno y coseno a incluir sería \(K=m/2\) donde \(m\) es la frecuencia de la estacionalidad de la serie.
Pronostico con regresión lineal
El pronostico del modelo estimado no es más que usar los \(\hat{\beta}\) encontrados por MCO sobre una nueva observación de \(X\)\[
\hat{y_t} = \hat\beta_{0} + \hat\beta_{1} x_{1,t} + \hat\beta_{2} x_{2,t} + \cdots + \hat\beta_{k} x_{k,t}
\]
Además un intervalo de confianza al 95% para le pronostico (asumiendo normalidad en los errores del modelo) se define como : \[
\hat{y} \pm 1.96 \hat{\sigma}_e\sqrt{1+\frac{1}{T}+\frac{(x-\bar{x})^2}{(T-1)s_x^2}}
\]
Pronostico ex-antes y ex-post
Pronostico ex-antes : Son aquellos que se realizan únicamente con los datos que están disponibles previamente de los predictores del modelo.
Pronostico ex-post : Son aquellos que se realizan utilizando información posterior o nueva sobre los predictores del modelo.
Otra forma de hacer pronosticos es con la construcción de escenarios, donde suponemos que las variables de nuestro modelo toman un rango de valores esperados en el futuro. Esto es útil para tomadores de desición si quireren ver valores esperados de la serie bajo supuestos y hacer comparaciones.
Ejercicio 5:
Hacer una regresión del consumo con el ingreso en el modelo macroeconómico de Estados Unidos. El ministro de economía del país está interesado en saber cúal será el consumo esperado en el siguiente trimestre en tres escenarios : si el ingreso se matienen en sus valores promedios historicos, se reduce 0.2 puntos porcentuales o aumenta 0.8 puntos porcentuales respecto al promedio historico.
\(\Delta C_t\) es el cambio porcentual en el gasto real de consumo personal.
\(Ingreso_t\) es el cambio porcentual en el ingreso personal real disponible.
Usar el archivo US_change.csv para modelar.
Rendimiento y comparación de modelos para la predicción
Cuando tenemos varios modelos con diferentes covariables pero con la misma variabla dependiente podemos usar estadísticos de selección de modelos como :
\(R^2\) ajustado : Permite comparar modelos teniendo en cuenta el número de covariables en el modelo. Entre más alto sea, mejor es el modelo. Al ser un \(R^2\) penalizado, este es inferior que el \(R^2\).
\[
\bar{R}^2 = 1-(1-R^2)\frac{T-1}{T-k-1}
\]
Donde \(T\) son el número de observaciones y \(k\) el número de covariables en el modelo.
Criterio de información de Akaike : Permite comparar modelos, donde el valor más pequeño es el mejor modelo.
Para hacer pronosticos, nos intereza verificar la presición del modelo usando los errores de estimación \(\varepsilon_t = y_t - \hat{y}_t\). Para ello, verificamos qué tanto se equica mi modelo con los datos, donde vamos a seleccionar el modelo con el menor error posible. En la práctica necesitamos dividir la serie en un segmento para estimar el modelo y otro para verificar el rendimiento.
Promedio absoluto de los errores escalados (MASE) : \[
q_{j} = \frac{\displaystyle e_{j}}{\displaystyle\frac{1}{T-1}\sum_{t=2}^T |y_{t}-y_{t-1}|} \\
\text{MASE} = \text{mean}(|q_{j}|)
\]
Raíz del error cuadrárico escalado (RMASE= : \(\text{RMSSE} = \sqrt{\text{mean}(q_{j}^2)}\)
Una técnica más avanzada es la validación cruzada (CV), donde el objetivo es mover el punto de cohorte del modelo e ir verificando el rendimiento en el conjunto de datos posterior. Luego se promedia un estadístico de precesición del pronostico (usualmente MAE):
Elimine la observación \(t\) del conjunto de datos y ajuste el modelo utilizando los datos restantes. Luego, calcule el error (\(e_{t}^*=y_{t}-\hat{y}_{t}\)) para la observación omitida. (Note que este no es el mismo valor que el residuo, ya que la t-ésima observación no se utilizó para estimar el valor de \(\hat{y}_{t}\)).
Repita el paso 1 para cada \(t=1,\dots,T\).
Calcule el MSE (Error Cuadrático Medio) a partir de \(e_{1}^*,\dots,e_{T}^*\). A este resultado lo llamaremos CV.
Ejercicio 6:
Usando el modelo macroeconómico del cambio en el consumo de Estados Unidos, compara el rendimiento de todas la combinaciones (16 modelos en total) de covariables disponibles en términos de \(Ingreso_t\), \(Industria_t\), \(Ahorro_t\), \(Desempleo_t\).
Calcular : AIC, BIC, CV, RMSE, MAE, MAPE, MASE para cada modelo.
Finalmente, compare los resultados y seleccione el mejor modelo.